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Abstract. Analytic solutions F{v,t) of the nonUnear Boltzmann equation in d- 
dimensions are studied for a new class of dissipative models, called inelastic repul- 
sive scatterers, interacting through pseudo-power law repulsions, characterized by 
a strength parameter v, and embedding inelastic hard spheres {f = 1) and inelas- 
tic Maxwell models {v — 0). The systems are either freely cooling without energy 
input or driven by thermostats, e.g. white noise, and approach stable nonequilib- 
rium steady states, or marginally stable homogeneous cooling states, where the 
data, VQ{t)F{v,t) plotted versus c — v/vo{t), collapse on a scaling or similarity 
solution /(c), where vo{t) is the r.m.s. velocity. The dissipative interactions gener- 
ate overpopulated high energy tails, described generically by stretched Gaussians, 
/(c) ~ exp[— /3c''] with < 6 < 2, where h = v with > in free cooling , and 
b = 1 + with v > when driven by white noise. Power law tails, /(c) ~ 1/c"^'', 
are only found in marginal cases, where the exponent a is the root of a transcen- 
dental equation. The stability threshold depend on the type of thermostat, and is 
for the case of free cooling located at = 0. 

Moreover we analyze an inelastic BGK-type kinetic equation with an energy de- 
pendent collision frequency coupled to a thermostat, that captures all qualitative 
properties of the velocity distribution function in Maxwell models, as predicted by 
the full nonlinear Boltzmann equation, but fails for harder interactions with v > 0. 

1 Introduction 

Classic kinetic theory [1-5] deals with elastic particles with energy conserving 
dynamics. The system is described by the single particle distribution func- 
tion, whose time evolution is governed by the nonlinear Boltzmann equation. 
The asymptotic states of such systems follow the universal laws of thermody- 
namics, and the distribution function is the Maxwell Boltzmann distribution. 
This scenario does not apply to dissipative systems, where energy is lost in 
inelastic interactions. 

In elastic systems the approach to asymptotic states is characterized by 
a kinetic stage of rapid relaxation in velocity space to a locally homogeneous 
equilibrium state, followed by a hydrodynamic stage of slow approach to a 
globally homogeneous equilibrium state. The time scale in the kinetic stage 
is the mean free time tmf between collisions. In the kinetic theory of inelastic 
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systems [6-12] the type of decay depends on the energy supply to the dissi- 

pativc system. Without energy supply there is first a kinetic stage of rapid 
relaxation on the time scale tmf to a locally homogeneous adiabatic state, 
the homogeneous cooling state, described by scaling or similarity solutions 
with a slowly changing parameter, at least for weakly inelastic systems. With 
energy supply the evolution is more similar to the elastic case with, however, 
equilibrium states replaced by non-equilibrium steady states. The velocity 
distributions in these adiabatic or steady states are very different from a 
Maxwell Boltzmann distribution. The subsequent stage of evolution involves 
transport phenomena and complex hydrodynamic phenomena of clustering 
and pattern formation [12,13]. 

The interest in granular matter in general has strongly stimulated new 
developments in the kinetic theory of granular fluids and gases, which show 
surprising new physics. A granular fluid is a collection of small or large macro- 
scopic particles, with short range repulsive hard core interactions, in which 
energy is lost in inelastic collisions, and the system cools when not driven. 
When rapidly driven, gravity can be neglected. The dynamics is based on 
binary collisions and ballistic motion between collisions, which conserve to- 
tal momentum. So these systems can be considered to be a granular fluid 
or gas. The prototypical model for these so-called rapid granular flows is a 
fluid or gas of perfectly smooth mono-disperse inelastic hard spheres, and its 
non-equilibrium behavior can be described by the nonlinear Boltzmann equa- 
tion [6-12]. The inelastic collisions are modeled by a coefficient of restitution 
a (0 < a < 1), where (l — a^) measures the degree of inelasticity. 

This review focuses on the first stage of evolution, and studies the velocity 
distribution F (v, t) in spatially homogeneous states of inelastic systems. For 
that reason most of the citations, given in this article, only refer to kinetic 
theory studies of F{v,t), and not to studies of transport properties. The 
revival [14-26] in kinetic theory of inelastic systems has been strongly stim- 
ulated by the increasing sophistication of experimental techniques [27,28], 
which make direct measurements of velocity distributions feasible in non- 
equilibrium steady states. In this review we also include inelastic general- 
izations [29] of the classical repulsive power law interactions [2 4], which 
embed both the inelastic hard spheres (i^ = 1), as well as the recently much 
studied [30-40] inelastic Maxwell models {ly = 0) in a single class of models, 
parametrized by an exponent v. This exponent characterizes the dependence 
of the collision frequencies on the energy of impact at collision. 

In fact, the kinetic theory for such models is of interest in its own right, 
as the majority of inter-particle interactions in macroscopic systems involve 
some effects of inelasticity. Our goal is to expose the generic and universal 
features of the velocity distributions in dissipative fluids, and to compare 
them with conservative fluids to highlight the differences. 

A classical problem in kinetic theory is the possibility of overpopulated 
high energy tails in velocity distributions [41,42], as many physical and chem- 
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ical processes only occur above a certain energy threshold. Consequently, this 

old problem has received a new stimulus through a large amount of recent 
theoretical and experimental studies on tail distributions in many particle 
systems with inelastic interactions. Prom the point of view of kinetic theory 
the intriguing question is, what is the generic feature causing overpopulatcd 
tails, possibly even power law tails, in systems of inelastic particles, how 
does the overpopulation depend on the scattering cross sections, and on the 
different forms of energy input. The generic feature is the mechanism for 
overpopulation, and not the specific shape of the tails. 

Finally, from the point of view of nonequilibrium steady states, the struc- 
ture of velocity distributions in clastic and dissipative systems, including the 
high energy tail, is a subject of continuing research, as the universality of the 
Gibbs' state of thermal equilibrium is lacking outside thermal equilibrium, 
and a possible classification of generic structures would be of great interest 
in many fields of non-equilibrium statistical mechanics. 

The plan of the paper is as follows: in Section 2 we discuss a simple inelas- 
tic BGK- or single-relaxation-timc model [43] to illustrate the phenomenon of 
power law tails. The exponent in the algebraic tail depends qualitatively in 
the same manner on the degree of inelasticity as in 2- and 3-dimensional 
Maxwell models. In Section 3 the nonlinear Boltzmann equation is con- 
structed for inelastic generalizations of the classical repulsive power law po- 
tentials, to which we refer as Inelastic Repulsive Scatterers or IRS-models. 
This is done for freely cooling as well as for systems driven by thermostats 
or heat sources. Section 4 gives a systematic analysis for the energy balance 
equation, it derives the nonlinear integral equation for the scaling or simi- 
larity solution, denoted by / (c) , and presents an asymptotic analysis of the 
high energy tails in the form of stretched Gaussians, / (c) ^ exp [—/3c''] where 
b < 2. The method used can only be applied to IRS-models with > 0, where 
the exponent b = b (v) is found as a simple function of ly. The case of freely 
cooling Maxwell models (i' = 0) forms a borderline case, discussed in Section 
5. Here algebraic tails, /(c) --^ l/c''"'"", are found, where the exponent a is 
determined by a transcendental equation. It yields a = a (a) as a function of 
the degree of inelasticity. We end with some perspectives and conclusions. 

2 Inelastic BGK Model 

2.1 Kinetic equations 

The goal of this section is to present in the nutshell of a simple Bhatnagar- 
Gross-Krook (BGK) model [43] a preview of many of the qualitative features 
of velocity relaxation in homogeneous systems. 

A crude scenario for the relaxation without energy input suggests that the 
system will cool down due to inelastic collisions, and the velocity distribution 
F {v, t) will approach a Dirac delta function 5 {v) as i ^ oo, while the width or 
r.m.s. velocity vq (t) of this distribution, defined as (f^) = ^dv^, is shrinking. 
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With a constant supply of energy, the system can reach a non-equihbrium 
steady state (NESS). 

To model this evolution we use a simple BGK-type kinetic equation, in- 
troduced in 1996 by Brey et al. [44] , 

dtF {v, t) - DVlF {v, t) - -oj, (t) [F (v, t) - Fo (v, t)] . (1) 

We have added a heating term, —DV^F {v, t) to the usual BGK equation 
which represents the heating by a white noise of strength D ^. Here the 
mean collision frequency, lij,^ (t) = 1/tmf, is a function of the r.m.s. velocity 
fo it), chosen as = Vq, in preparation of section 3.1. The kinetic equation 
describes the relaxation of F {v, t) with a time-dependent rate {t) towards 
a Maxwellian with a width proportional to vq (t), defined by 

Fo {v, t) = (Vttoto) oxp - {v/av^f = {av^y^ (j) (c/a) , (2) 

where c = v/vq- The constant a (0 < a < 1) is related to the inelasticity, 
7 = I (l — a^), of the model, and the totally inelastic limit (a — > 0) is 
ill-defined in this model, as the mean energy is divergent for a = 0. 



2.2 Free Cooling (D = 0) 

The cooling law of the mean square velocity (w^) = ■^dv^, or the granular 
temperature T = VQ,is obtained by applying J dvv^ (...) to (1) with lOi, = Vq. 
The result is dfVo = — T^o^^, yielding 

^^o(^)=^^o(0)/[l + ^/7^^'o(0)]^/^ (3) 

The result is a homogeneous cooling law, T t~'^^'^, which agrees with 
Haff's law [9] for i/ = 1, corresponding to inelastic hard spheres. Note that 
for negative u the homogeneous cooling law takes the form, 

vo{t) = vo{0)[l-Wht/vl''^{0)]'^^'^\ (4) 

i.e. for t > tg = Vq'^/WIj the r.m.s. velocity and the mean kinetic energy 
become negative, which is unphysical, and so are the BGK-models with < 0. 

As indicated in the introduction, the long time behavior of F {v, t) in free 
cooling is determined by a scaling or similarity solution of the form F (v,t) = 
Vq"^ {t) f {v I Vq {t)) . We insert this ansatz in (1), eliminate dtVQ = —7^0^^, 
and obtain the scaling equation, 

^ For a systematic discussion of driven systems, see Section 3.3 
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where a is defined as, 

a = 1/7 = 2/ (1 - . (6) 

Wc also note that the scaling equation is independent of u. The exact solution 
of this equation is: 

/ W = ^ + ^ (^) [duu^^^-^e^ [-u^/a-] . (7) 
The integration constant A is fixed by the normalizations, 

j dc{l,c^}f{c) = {l,\d}, (8) 

where dc = ndc'^~^dc with 12^ = / F {d/2) being the surface area of a 
rf— dimensional hyper-sphere. The normalization integral converges only near 
c ~ if ^ = 0. Then the solution (7) with A = is identical to the scaling 
form, obtained in [44] for j/ = 1. 

For velocities far above thermal, i.e. 1, the solution has a power law 
tail, 

with exponent a = I/7 = 2/ (l — a^), such that (c^) is bounded for a > 0. 
The exact solution (7), including its high energy tail, is independent of the 
exponent i/, that determines the energy dependence of the mean coUision 
frequency lu,^ ~ Vq in the BGK model. 

As we shall see in Section 5, a similar heavily overpopulated tail, / (c) ^ 
l^(.d+a ^[^ly d > 1, will also be found in freely cooling Maxwell model with 
Wo = 1 (i^ = 0). There the exponent a (a) takes in the elastic limit (a 1) 
the form a ~ I/70 = 4c?/ (l — a^). However, in the general class of IRS- 
models with collision frequency loi^ ^ Uq with i> > the tails are not given 
by power laws, but by stretched Gaussians, f (c) ^ cxp [—/3c''] with < 6 = 
b (f) < 2. These models will be introduced and discussed in Section 3.2. 



2.3 NESS {D ^ 0) 

Next we extend the result of [44] to an inelastic BGK-model with ly > 0, 
driven by white noise. By applying J dvv"^ to (1) we obtain the temperature 
balance equation at stationarity, 

dtvl =4D- 27o<+2 = 0, (10) 

where the collisional dissipation, 270^0^^, is compensated by energy input 
from the external white noise. Also note that for u < —2 the fixed point 
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solution vo (oo) in (10) still exists, but it is unstable [29]. If vo (0) < vq (oo), 

then vq (t) vanishes a.s t ^ oo, and if vq (0) > vq (oo) then vq (t) diverges ^. 

To obtain the solution F {v, oo) of (1) in the NESS we rescale F {v, oo) = 
Vq''' (oo) / {v/vo (oo)) to the standard width (c^) = ^d, substitute the rescaled 
form in the kinetic equation, and eliminate D, using the stationarity condition 
(10) as well as (6). This yields the rescaled equation in universal form, 

1 d J T d „ , , ^ „ , . 2a 



c<^-i dc 



'|/(c)-2a/(o) = --*(£). (II) 



where the normalizations (8) are imposed. The O.D.E. shows that /(c) is 
independent of the noise strength, D, and does not contain any dependence 
on the initial data. This equation can be solved exactly, and more details will 
be published in [29]. However, for the purpose of this section, we only want 
to extract from the differential equation (11) the asymptotic form of / (c). In 

that case we may neglect in (11) the inhomogcncity (c/a) ^ cxp — {c/af 

for c » a, and find the asymptotic solution for the BGK-model driven by 
white noise, i.e. 

/(c)^exp[-/3c^] 



& = 1; (3 = V2^ = 2/Vl-a2. 

The constant (3 is independent of the parameter v. The exponentially decay- 
ing high energy tail is also a 'stretched' Gaussian, which is overpopulated 
when compared to a Maxwellian, but the overpopulation is much less heavy 
than is the freely cooling case (9) with an algebraic tail. 

As wc shall see in Section 4.3, a similar exponential high energy tail will 
be found in the white noise driven Maxwell model {v = {): 6=1; (3 = 
■\/8/ (1 — a^)), but not in the general class of IRS-models, where the stretch- 
ing exponent b takes a value in the interval < 6 = 6 [u) < 2. 



3 Basics of inelastic scattering models 

3.1 Boltzmann equation as a stochastic process 

The nonlinear Boltzmann equation for dissipative interactions in the homo- 
geneous cooling state can be put in a broader perspective, that covers both 
elastic and inelastic collisions, as well as interactions where the scattering 
of particles is described either by conservative (deterministic) forces, or by 
stochastic ones. To do so it is convenient to interpret the Boltzmann equation 
as a stochastic process, similar to the presentations in the classical articles 
of Waldmann[l], and Uhlenbeck and Ford [5], or in Ulam's stochastic model 

^ Thanks are due to E. Trizac and A. Barrat for pointing out to us that the energy 
balance equations in this article have stability thresholds that are z/— dependent. 
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[45] showing the basics of the approach of a one-dimensional gas of elastic 

particles towards a MaxwcUian distribution. 

Consider a spatially homogeneous fluid of elastic or inelastic particles 
in d— dimensions, specified by their velocities {v,w,---), and interacting 
through binary collisions, {v,w) {v',w'), that arc described by transi- 
tion probabilities. To describe fluids out of equilibrium the total momentum 
G = ^{v + w) needs to be conserved in a binary collision. The outgoing or 
post-collision velocities (v' , w') can be parametrized in terms of the incoming 
velocities {v, w), and an impact (unit) vector n, that is chosen on the surface 
of a unit sphere with a certain probability, proportional to the collision fre- 
quency A {g, n) = A (g, g • n), that in general depends on the relative speed 
g = \v ~ w\ of the colliding particles, and the angle between g and n, where 
d denotes a unit vector. 

In this article we consider the simplest case of modeling the inelastic colli- 
sions through a velocity- independent coefficient of restitution a (0 < a < 1), 
where the component g|| = g • n is not only reflected, as in elastic collisions, 
but also reduced in size by a factor a, i.e. 

g||' = -ag||. (13) 

The components g± = g — g\\n, orthogonal to n, remain unchanged. More 

explicitly, with the help of momentum conservation we find for the post- 
collision velocities resulting from the direct collisions {v,w) {v',w') = 
{v*,w*), 



V* =v — ^{1 + a) g ■ nn 
w* =w + ^{1 + a)g ■ nn. 



(14) 



The corresponding energy loss in such a collision is, 

AE=^ {v*^ + w*^ -v^- w^) = -Hi - (9 ■ , (15) 

where (l — a^) measures the degree of inelasticity. The value a = 1 describes 
elastic collisions. In the special case of hard spheres the impact vector n is 
the unit vector along the line of centers of the colliding spheres at contact, 
as illustrated in Figure 1. 

The system above can be described by an isotropic velocity distribution, 
F{v,t) = F{\v\,t), as long as one only considers isotropic initial distribu- 
tions. Its time evolution is given by the Boltzmann equation, 

dtF{v,t) = Iiv\F) = j dwdv'dw' j dn[W {v,w\v' ,w':n) 

xF{v',t)F{w',t) - W {v' ,w'\v,w;n) F {v,t) F {w,t)] . (16) 



Here W {v' ,w'\v,w;n) is the transition probability per unit time that the 
incoming pair state {v,w) at impact vector n is scattered into the outgoing 
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Fig. 1. Elastic (a = 1) and inelastic (a < 1) scattering of hard spheres, where 

the parallel component gy is reflected as in (13). In elastic collisions the scattering 
angle x = cos~^ io ■ o') = ^ — with (j) = cos^^ {g ■ n) and in inelastic ones 
X = cos~^ {g ■ g*) = ir — cj) — (f>* with 0* = cos~^ {g* ■ n). 

pair state (v' , w'). The loss term is the sum over all parameters of the direct 
collisions {v,w) {v\w') at fixed v, including a sum over impact vectors 
n. Similarly the gain term is the sum over all parameters of the restituting 
collisions, (f', w') {v, w) at fixed v. The transition probability for the scat- 
tering event {v,w) {v',w'), obeying the inelastic reflection law (13) and 
momentum conservation, is in general proportional to the collision frequency 
A{g,n), and contains delta functions, selecting the allowed collisions, 

W {v', w'\v, w;n)=A {g, n) 5^ {G' - G) 5^''-'^ {g'^ - g^) 5 (g,,' + ag,,) . 

(17) 

The faster the collision frequency increases at large impact velocities, the 
more rapidly the high energy tail of F {v, t) relaxes relative to the bulk values 
of F [v, t) with V in the thermal range. 

For conservative interactions the total energy is conserved, as well as 
total momentum and total number of particles. Moreover (17) shows that 
the transition probability is symmetric for a = 1, 

W {v', w'\v, w;n) = W {v, w\v', w'; n) . (18) 

This means that the transition probabilities for elastic binary collisions sat- 
isfy the condition of detailed balance. Energy conservation and the detailed 
balance relation in combination with the if— theorem guarantee that the 
Maxwellian velocity distribution is approached at large times. 
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Once the detailed balance relation is obeyed, it is trivial to prove the 
-ff— theorem. Defining the if— function or entropy S (t) = —H {€) as, 

H{t)= j dvF{v,t)lnF{v,t), (19) 

wc obtain its decay rate by applying J dvhiF {v,t) to (16), then symmetriz- 
ing over V ^ w and v' w' , and subsequently over {v,w) ^ {v' ,w'). The 
result is, 

dtH it) = jj dvdwdv'dw' J dnW [Y - X] In (X/Y) < 0, (20) 

where X = F (v) F {w) and Y = F {v') F {w') and the inequality follows from 
[Y - X] In (X/Y) < 0. The equality sign holds if and only iiX = Y. This 
implies that H decreases monotonically and becomes stationary only when 
F approaches the Maxwellian. 

The reason for reviewing these 'obvious' properties in elastic systems is 
that several of them, such as energy conservation, the detailed balance rela- 
tion, the iJ— theorem, and the approach to a stationary Maxwellian velocity 
distribution no longer hold in dissipative systems (a < 1), as we shall see. 



3.2 Boltzmann equation in standard form 

The standard form of the collision term I {v\F) in the Boltzmann equation for 
elastic interactions (a = 1) contains the differential scattering cross-section 
a {g, x). It can be calculated from the pair potential V (r) [4] , and depends on 
the relative speed g — \g\ and the scattering angle x = — 2ip, where x ^nd 
(f are defined in Figure 1. The collision frequency is then given by A (g, n) = 
g(j{g,x)- For repulsive power law potentials V (r) ~ r~* the collision fre- 
quency is A {g, n) = g" A {g ■ n) with i' = 1 — 2 {d — 1) /s in d— dimensions [46]. 
By definition Maxwell molecules have a collision frequency, which is indepen- 
dent of g, corresponding to = or s = 2 (d — 1). For hard spheres, s ^ oo or 
p = 1 with A (g, n) = \g ■ n\ a'^~^ = |g|| | a'^~^, and a is the hard sphere diam- 
eter. Using these relations the collision term in (16) can be reduced to the clas- 
sical Boltzmann equation for elastic energy-conserving collisions. For inelastic 
hard spheres the collision rate is again given by A (g, n) = |g|| | cr'^^^. In gen- 
eral, a given positive function A (g, n) defines a stochastic scattering model, 
and in particular the choice, A{g,n) — A{g ■ n), defines inelastic Maxwell 
models [24-26]. As a simple realization of a dissipative scattering model we 
consider the IRS-models with a collision frequency A (g, n) = Aq |g • n\'^ with 
1/ >0. This class includes inelastic hard spheres (i/ = 1) and inelastic Maxwell 
models = 0). In the remainder of this article we restrict ourselves to this 
class of models. 

To reduce the Boltzmann equation to its standard form we use the tran- 
sition probabilities for dissipative interactions (17) with A (g,n) = Aq |g|| | . 
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Consider first the loss term in (16), insert the relation dv'dw' = dg'j_dg\\' dG' , 
and carry out the integrations over the delta functions. The result is, 

{dtF (^;))ioss = -Jdwl dnAo |g|| j'' F (w) F (v) . (21) 

The gain term involves the transition probability for the restituting collisions 
{v', w') {v, w), obtained from (17) by interchanging primed and unprimed 
velocities, 

(atF(^;))gain = j dwdg'^dg\{ dG' j dnAo IguT^''^' (G - G') 

X S^"-'^ [g^ - g'J (^g||' + ig||^ F {w') F {v') (22) 

= J dw j dn {1/ a) Ao \gii / a]"" F{w**)F{v**). 

In the second integral we have carried out the integrations over the primed 
velocities, and used the following relations for the restituting velocities, 

v' =G' + + igii'n = G+\g^- ^gpn 

=v-\{^ + 7i)9-nn = v** (23) 

w' =w + I (l + ^) 9 • nn = w**. 

In the first equality v' has been expressed in center of mass and relative veloc- 
ities. In the second equality we have used the inelastic collision law (13) and 
conservation of total momentum, and the very last equality defines the resti- 
tuting velocities, {v**,w**). They are the incoming velocities that result in 
the scattered velocities {v, w), described by the inverse of the transformation 
(14). 

The space-homogeneous Boltzmann eqiiation in its standard form is then 
obtained by combining (21) and (22) with A (g, n) = \g ■ n\'^ to yield, 



dtFiv) =I{v\F) 

I{v\F) =jjdw Ig-nl" 



^ F (v**) F (w**) - F (v) F {w) 



(24) 



where /„(••) = {1/Od) / dn{- ■■) is an average over a rf— dimensional unit 
sphere, and we have absorbed constant factors in the time scale. Here u = \ 
corresponds to inelastic hard spheres and = to inelastic Maxwell models. 
Velocities and time have been dimensionalized in terms of the width and the 
mean free time of the initial distribution. Moreover, the Boltzmann equation 
obeys conservation of particle number and total momentum, but the aver- 
age kinetic energy or granular temperature, T ~ (^'^), decreases in time on 
account of the dissipative collisions, i.e. 

j dv (1, V, v^) F {v, t) = (1, 0, \dvl it)) , (25) 
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where vq (t) is the r.m.s. velocity. The inelastic scattering models with colli- 
sion frequency A ^ (i^ > 0), arc the inelastic analogs of the deterministic 
models with repulsive power law potentials, V (r) ~ with 2{d — 1) < 
s < oo, introduced at the start of this section. The inelastic case with v = 2 
corresponds to an exactly solvable stochastic scattering model with energy 
conservation, known as Very Hard Particle model [46]. 

The most basic and most frequently used model for dissipative systems 
with short range hard core repulsion is the Enskog-Boltzmann equation for 
inelastic hard spheres in d— dimensions [6], which simplifies in the spatially 
homogeneous case to (24) with v = 1. Recently, inelastic Maxwell models have 
been studied extensively. Ben-Nairn and Krapivky [24] introduced the one- 
dimensional version of (24) with [d= = 0), and Bobylev et al [25] have 
introduced a three-dimensional Maxwell model with A {g, n) = Ao\g ■ n\. The 
Maxwell models in (24) for general d were first considered in [32,34]. 

3.3 Cooling and driven systems 

Homogeneous cooling and scaling: 

An inelastic fluid without energy input will cool down due to the coUisional 

dissipation in (15). In experimental studies of granular fluids energy has to be 
supplied at a constant rate to keep the system in a non-equilibrium steady 
state, while in analytic, numerical and simulation work freely cooling sys- 
tems can be studied directly. Without energy input the velocity distribution 
F {v,t) will approach a Dirac delta function S{v) ast ^ oo, and all moments 
approach zero, including the width vq (t). 

However, an interesting structure is revealed when velocities, c = v/vq (t), 
are measured in units of the instantaneous width vq {t) , and the long time 
limit is taken while keeping c constant, the so-called scaling limit. Monte 
Carlo simulations [30] of the Boltzmann equation suggest that in this limit 
the rescaled velocity distribution of the homogeneous cooling state can be 
collapsed on a scaling form or similarity solution /(c). These observations 
seem to indicate that the long time behavior of F (v, t) in freely cooling sys- 
tems approaches a simple, and to some extent universal, form / (c), which is 
the same for different initial distributions. Such scaling or similarity solutions 
have the structure. 



where c = v/vq (i) is the scaling argument. Then / (c) satisfies the normal- 
izations. 



F{v,t) = {vo{t))-''f{v/vo{t)), 



(26) 




(27) 
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on account of (25). Substitution of the scaling ansatz (26) in the Boltzmann 
equation (24) then leads to a separation of variables, i.e. 



where 7 is a separation constant. 

The question is then: Can one determine physically acceptable solutions? 
Here a short history. The first kinetic model with dissipative interactions, 
that has been solved exactly for the case of free cooling, and exhibits a heavily 
overpopulated high energy tail, /(c) ^ l/c''+", is the inelastic BGK model, 
discussed in Section 2. Asymptotic solutions (c ^ 1) of the scaling equation 
for inelastic hard spheres [16,17] have predicted the existence of exponential 
high energy tails, f (c) ~ exp [—/3c]. We talk about tails, over-populations or 
stretched Gaussians in / (c) , when the ratio of / (c) and a Gaussian is an 
increasing function of c at c ;:g> 1. The predictions about high energy tails 
were later confirmed in great detail by Monte Carlo simulations of the long 
time solutions of the nonlinear Boltzmann equation for inelastic hard spheres 
[20,21], as well as in further analytic work [18]. 

Similarity solutions for freely cooling inelastic Maxwell models [24,25] 
where first studied in terms of the scaling variable, (l — a^) t, relating large 
times and small inelasticities [25], but the solutions obtained turned out to 
be unphysical, i.e. non-positive. The first exact positive similarity solution 
was found by Baldassarri et al. [30] for the one-dimensional Maxwell model. 
It shows a surprisingly strong high energy tail of algebraic type, ~ 1/c"'. Us- 
ing Monte Carlo simulations of the Boltzmann equation, these authors also 
showed that rather general initial distributions approach towards this scaling 
form for long times. Algebraic tails, /(c) ^ l/c'*+", for Maxwell models in 
higher dimensions have been obtained analytically in [32 36], and the actual 
approach in time towards such scaling forms has been studied analytically in 
[36,37] for general classes of initial distributions. 

Driving and steady states: 

Heating may be described by applying an external stochastic force to the 
particles in the system, or by connecting the system to a thermostat, which 
may be modeled by a frictional force. For example, the friction force 'yv - here 
with a negative friction coefficient (—7) - is called a Gaussian thermostat. 
Complex fliuds (e.g. granular) subject to such forces can be described by the 
microscopic equations of motion for the particles, r, = Vi, and Vi — Ui -\- ^^ 
{i = 1, 2, • • • ), where a, and |j are respectively the systematic and random 
forces per unit mass. Here at contains frictional forces, which may depend 
on velocity, and in the present case represents external noise (modeling 
energy input), which is taken to be Gaussian white noise with zero mean, 



I (c|/) =7 {df + cdf/dc) = 7Ve • cf 

Vo=- JVq^ , 



(28) 
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and variance, 



^i,a it) ^j^fs (to = 2DSijSa,i36 [t - t') , (29) 

where a, (3 denote Cartesian components, and D is the noise strength. The 
Boltzmann equation for system driven in this manner takes the form, 

dtF{v)+TF{v) = I{v\F) 



(30) 



A detailed derivation on how to include frictional and stochastic forces in 
kinetic equations can be found, for instance, in [47,14,17]. 

When energy is supplied at a constant rate, driven systems can reach a 
NESS. Again there is the question about universality of these asymptotic 
states as t — > oo. Does the NESS depend on the inelasticity, on the type of 
thermostat, and on the initial distribution? The basic idea to show universal- 
ity is always essentially the same. Rescale the velocity distribution F {v, oo) 
by measuring velocities in terms of their typical size, i.e. the width vq (oo), 

F(v,^) = (vo{^)r''f{v/vo{^)), (31) 

and analyze the scaling form / (c). The scaling equation for the NESS func- 
tion / (c) will be analyzed in Sections 4.2 and 4.3. Here we make some com- 
ments about what is known. The scaling solutions for the NESS show again 
overpopulated tails in the form of stretched Gaussians, albeit with a larger 
stretching exponent b than in free cooling. The scaling solution for inelas- 
tic hard spheres driven by white noise was first analyzed in [17], predicting 
/(c) ~ exp [—/3c''] with b = 3/2. The experiments of Rouyer and Menon 
[28,27], and of Aranson and Olafsen [28] seem to confirm the stretched Gaus- 
sian behavior with b = 3/2. Maxwell models, driven by white noise, were first 
studied in [24,26], but did not give any predictions about overpopulated high 
energy tails. The high energy tails for Maxwell models, driven by white noise, 
were predicted in [35] to have exponential tails, / (c) ^ exp [—/3c]. A number 
of more detailed analytical and numerical studies about this driven Maxwell 
model in one dimension [38-40,48] have appeared as well. Moreover, as will 
be shown in Section 4.2, the integral equation for the scaling form / (c) in 
free cooling is identical to the integral equation for the NESS distribution for 
a special thermostat, provided F {v, oc) is also rescaled to the same width as 
in (27). How to determine the scaling form /(c) in free cooling and NESS 
will be described in subsequent sections. 



3.4 Qualitative analysis 

In order to illustrate the rich behavior of the inelastic systems, we start with 
the Boltzmann equation for the one-dimensional inelastic Maxwell model (24) 
{d= 1,^ = 0), where F (v, t) satisfies. 



dtF {v) - D^lF {v) = I {v\F) , 



(32) 
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and the coUision term has the form, 

I{v\F) =jdw 

= -F(v) + - [ duF (u) F I — 
PJ VP 



^F (v**) F (w**) - F {v) F {w) 
1 f . f V — qu 



(33) 



All velocity integrations extend over the interval (— oc', +oo). Here the outgo- 
ing velocities {v,w), and the incoming ones {v**,w**) are according to (23) 
related by, 

V = qv + pw ; w = pv + qw (34) 

with p = 1 — q = i (1 + a). By changing integration variables w v** = u 
with dw = {a/p) du, and using the relation w** = {v — qu) jp one obtains the 
second equality in (33). The normalization of mass and mean square velocity 
(^^) ~ 5^0 given by (25) for d=\. The temperature balance equation is 
obtained from (32) as, 

dtvl = 4D- 2pqvl (35) 

and describes the approach to the non-equilibrium steady state (NESS) with 
a width Vq (oo) = 2D/pq, where the heating through random forces, ~ D, is 
balanced by the collisional losses. 

To understand the physical processes involved we first discuss in a qual- 
itative way the relevant limiting cases. Without the heating term (£) = 0), 
equation (32) reduces to the freely cooling inelastic Maxwell model. 

If one takes in addition the elastic limit (a — > 1 or g — > 0), the collision 
laws reduce in the one- dimensional case to v** = w, w** = v, i.e. an exchange 
of particle labels; the collision term vanishes identically; every F {v, t) = F (v) 
is a solution; there is no randomization or relaxation of the velocity distri- 
bution through collisions, and the model becomes trivial at the Boltzmann 
level of description, whereas the distribution function in the presence of in- 
fijiitesirnal dissipation (a ^ 1) approaches a Maxwellian. 

If we turn on the noise {D ^ 0) at vanishing dissipation (g = 0), the col- 
lision term in (32) vanishes, and the granular temperature follows from (35) 
as Uq {t) = Vq (0) -|- ADt, and increases linearly with time. With stochastic 
heating and dissipation (even in infinitesimal amounts) the system reaches a 
NESS. 



3.5 Comments 

Spatial dependence: 

When spatial dependence of the distribution function F (r,v;t) is relevant, 
the collision term in the Enskog-Boltzmann equation must be slightly mod- 
ified, namely the angular integration over the full solid angle should be 
replaced by 2 J^6 {—g ■ n), where the unit step function 9 (x) restricts the 
n— integration to the pre-coUision hemisphere with g ■ n < (see Figure 1). 
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In the spatially uniform case both representations are identical because the 
restituting velocities (23) are even functions of n. 

Origin of Maxwell models: 

In several papers [25,32,36,38] inelastic Maxwell models have been introduced 
as more or less ad hoc mathematical simplifications of the nonlinear collision 
term for inelastic hard spheres. This has been done by replacing the relative 
velocity g in the hard sphere collision frequency A{g,n) = \g ■ n\ by its mean 
value, (g) A (n) ~ vq {t) A (n), where ?;o (t) is the root mean square velocity. 

This proc;edure guarantees that the homogeneous cooling law for inelastic 
Maxwell model constructed in this way, is identical to the one for inelastic 
hard spheres, and given by Haff's law [9], T (t) ~ t~'^. The construction of 
inelastic Maxwell and IRS models, followed in the present article, is more 
in the spirit of [24], i.e. by defining the collision term through transition 
probabilities for the scattering process {v,w) ^ {v',w') with the proper 
constraints. 

Particles with pseudo-power law repulsion: 

A system of N inelastic hard spheres is a model with microscopic particles, 
characterized by positions and velocities, and interacting via well-defined 
force laws, that can be studied by means of MD simulations. However micro- 
scopic particles with dissipative inter-particle forces, such as inelastic Maxwell 
molecules and IRS-models, are stochastic models, only defined in AT— particle 
velocity space. Molecular Dynamics simulations can not be performed for such 
models. 

Violation of H— theorem: 

As the detailed balance symmetry and the iJ— theorem are lacking in dissipa- 
tive interaction models, there is no guarantee that the entropy S (t) = —H {t) 
is non- decreasing. In fact the solutions of the Boltzmann equation for such 
dissipative interaction models approach [36] for long times to a scaling form, 
defined in (26). By inserting such solutions into (19) and anticipating the 
cooling law (3) in the next section, one easily verifies that the entropy in the 
scaling state keeps decreasing as t becomes large, i.e. 

S{t) = -H{t)r^ - {d/v) In t. (36) 

This result is typical for pattern forming mechanisms in configuration space, 
where spatial order or correlations are building up, as well as in chaos theory, 
where the rate of irreversible entropy production is negative on an attractor 
[49]. Moreover, there is no fundamental objection against decreasing entropies 
in an open system in contact with a reservoir, which is here the energy sink 
formed by the dissipative collisions. The dynamics in A''— particle velocity 
space corresponds to a contracting flow dv*dw* = advdw where a < 1, 
where the probability is contracting onto an attractor. This is a well known 
phenomenon in chaos theory. 
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Modeling dissipation: 

There are many ways to model inelastic collisions that dissipate the relative 
kinetic energy of coUiding particles; e.g. by using Hertz' contact law [50,51], 
visco-elastic media [52], or coefficients of normal and tangential restitution 
[6,11]. The coefficients of restitution may also depend on the relative speed 
of the colliding particles [50] . 

For the IRS-models, studied in this article, the scattering laws (scattering 
angle, collisional energy loss) are independent of ^, i.e. they are the same for 
inelastic hard spheres, inelastic Maxwell models and for general IRS-models. 
Only the collision frequency A(g, n) of the IRS-models has the same energy 
dependence as the collision frequency of elastic particles interacting through 
repulsive power law potentials. 

4 Analysis of inelastic scattering models 

4.1 Homogeneous cooling laws 

We start with the freely evolving case and compute the cooling rate, defined 
through, dfVQ — —C,^ {t)vQ. This can be done by applying (/ dvv"^) to the 
Boltzmann equation (24), changing integration variables {v, w) {v**,w**) 
in the gain term, and using the relations dvdw = adv** dw** together with 
the energy loss AE per collision in (15). In general the above cooling equation 
is a formal identity, and does not provide a closed equation for Vo{t). However, 
if F{v, t) rapidly relaxes to a scaling form, then the subsequent time evolution 
of Vq — {2/d){v'^) is described by a closed equation for vo{t), as we will see, 
i.e. 



dtvl = I / dvv^I {v\F) = 11^1 dvdw \g ■ nl" AEF {v, t) F (w, t) 




(37) 



from which the cooling rate can be identified as, 



it) = -dtvl/vl = 270K.+2 vl (t) . (38) 



Here the coefficient 7o = 3^ (l — Oi^) measures the inelasticity. Moreover, the 
constant Ki, is defined as, 



Kv= j / dcdci\{c- ci) ■ n\'' f {c) f {ci) = f3y dcdci\c - ci}" f (c) f (ci) , 



Jn J 



(39) 



and the {d — 1) —dimensional angular integral is evaluated as. 
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In one- dimensional systems = 1 for all v, and in higher dimensions for 

the i^- values of interest here /J^ is convergent (z/ > 0). As a consequence of 
the scaling ansatz the coefficient time independent constant, that 

does depend on the unknown scaling form. For Maxwell models [u = 0) the 
cooling rate Co can be calculated explicitly from (40), and the constant K2 is 

^2 = i^/d) cif^ = 1 on account of (27). The result is, 

Co = -dtvl/vl = 270. (41) 

The r.m.s. velocity decays then as ?;o (t) = vq (0) exp [— 7ot]- 

The most important time scale in kinetic theory is the mean free time 
im/) which equals the inverse of the mean collision rate uji, {t), defined as the 
average of the collision frequency A {g,n) = \ {c — ci) ■ n\'^ over the velocities 
of the colliding pair in the scaling state, i.e. 

(i) = - y /loss (^1^^) = «-^o it) ■ (42) 
Comparison of (38) and (42) shows that both frequencies are related as, 

Cu {t) = 270 ^ (i) = 27,0;, (i) , (43) 

where 71/ Ki/ = 70K1/+2 are time independent constants. Using (42) and (43) 
the equation for the r.m.s. velocity becomes, dtVo = —^u^vVq'^^, yielding a 
solution, identical to (3) with 7 replaced by 7^. Consequently the granular 
temperature at large times decays asT = Vq t~'^/^ . For v = 1 one recovers 
the well known law of Haff [9] , describing the long time decay of the granular 
temperature in the homogeneous cooling state of inelastic hard spheres. A 
homogeneous cooling law similar to (3) with 1/ = 6/5 has been derived in [50], 
not on the basis of the Boltzmann equation (24) with an energy dependent 
collision rate cx g^^^ , but by using the Boltzmann equation for inelastic hard 
spheres {v = 1) with an energy dependent coefficient of restitution a (g) = 
1 — Ag"^/^. For Maxwell models {u — > 0) the r.m.s. velocity in (3) reduces to 
Vq (t) = Vq (0) exp [—70^], in agreement with (41). 

The homogeneous cooling law for the general class of IRS-models dis- 
cussed here, can be cast into a universal form by changing to a new time 
variable, the collision counter or internal tmic r of a particle, which represents 
the total number of collisions that a particle has suffered in the (external) 
time t. It is defined through the mean collision frequency, 

dr = (t) dt, (44) 

where C0v{t) ~ ^^{t). After inserting (3) in this equation the differential equa- 
tion can be solved to yield. 



VJiyT = In [H- V-fiyU!iy (O) t] , 



(45) 
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valid for all IRS-models. Combination of (45) and (3) shows that the r.m.s. 
velocity follows the universal homogeneous cooling law, 

^^0 (t) = vo (0) exp [-7^r] , (46) 

where 7^ is according to (43) and (39) proportional to the inelasticity 70, and 
depends on u through the collision integrals Ki,. We further observe that the 
relations (45) and (46) also holds for the inelastic BGK-model in Section 2 
with 7,y replaced by 7 = i (l — a^) in (6). 

So far we have been dealing with freely cooling systems. Next we address 
the balance equation for the granular temperature in driven cases, where the 
external input of energy counterbalances the coUisional cooling, and may lead 
to a NESS. We proceed in the same manner as for the free case, and apply 
(/ dw"^) to the Boltzmann equation in (30) with the result, 

dt (v^) = (t) (v^) +2{va)+ 2dD, (47) 

where the first term, J dvv"^! {v\F) = —(i, {t) {v"^), is obtained from (37), (38) 
and (25). The next two terms are obtained from the driving term in (30), i.e. 
J dvv^TF (v), by performing partial integrations. 

The most common ways of driving dissipative fluids [17,21,24,26] is by 
Gaussian white noise (WN) (a = 0;Z) 7^ 0), or by a Gaussian thermostat 
(GT) (a = "yv; D = 0), yielding for the balance equations. 



(27-C.(i))^o'(i) (GT) 
4D-C.(0«o(i) (WN). 



5t-o(i) = r;.' (48) 



Here the coUisional loss, — Ci/Wq, is counterbalanced by the heat, 27t;o, gener- 
ated by the negative friction of the Gaussian thermostat, or by the heat, 4D, 
generated by randomly kicking the particles. 

As t ^ oo the granular temperature will reach a NESS with T'(oo) = 
Vq (oo), where the r.m.s. velocity follows from (48) using (38), 

f(7/7o«.+2)'/'^ (GT) 
^°^°^^ = l(2D/7o«...)^/^-^) (WN). ^''^ 

With reference to the discussion in Section 2.3 we note that the GT-fixed 
point solution vq (oo) is attracting for u > 0, leading to a stable NESS, and 
unstable ior v < with vq (t) vanishing at t ^ oo if (0) < vq (oo), and 
diverging if ijq (O) > fo (oo). The case v = is m,arginally stable. For the 
WN-fixed point similar observations apply with the stability threshold i/ = 
replaced by v = —2. 

4.2 Scaling and non-equilibrium steady states 

Free cases: 

To investigate the existence of scaling solutions of the Boltzmann equation 
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for freely cooling inelastic systems, we substitute the scaling ansatz (26) into 
(24), to obtain an integral equation for the scaling form / (c). With the help 
of (44), (46) and (28) the left hand side of (24) becomes, 

dT\ ( ,„ , , d 



l.L.s. = -vov^^-^ ^-j|rf/(c) + c-/(c)^ 

= 7,0;, (r) (t) • cf (c) = 7o«.+2Vc • cf (c) ^"-^ (r) . (50) 

In the last equality we have used the relation j^k^ = implied by 

(43). The resulting integral equation for / (c) becomes, 

I (c|/) = ^0K„+2Vc ■ cf = ^i^c ■ cf. (51) 

With the help of (37) and (38) the second moment of the collision term, 072, 
can be expressed as, 

W2 = - j dcc^I (c|/) = rf7oK,+2. (52) 

Driven cases: 

We consider equation (30) for the NESS distribution in its rescaled form (31) 
with fg = Vq (00), driven by a Gaussian thermostat (GT), {a = "fv; D — 0} 
or by white noise (WN), {a = 0;£'^0}. The scaling equations for / (c) take 
the form, 

^Vc • cf = ]-m2Vc ■ cf (GT) 

1 (53) 

-72 _ -'^ „ V72 



^(ci/)=<^ 



= ~7r.^2Vif (WN). 



<+2 2d 



The first equality in (GT) and (WN) suggests that / (c) depends explicitly on 
7 or D. This is however not the case. The stationarity relation (49) combined 
with (52) shows in fact that the following expressions, 

072 = d-^v^" = 2dDvg'-''^^^ = rf7o'«!/+2, (54) 

are independent of 7 or D. So, we have used (49) to eliminate 7 and D, and 
to put it in the universal form, containing W2- 



4.3 Comments 

Equivalence free cooling - Gaussian thermostat: 

Comparison of integral equation (51) for free cooling and (53) for the Gaus- 
sian thermostat shows that both equations are identical, as first observed by 
Montanero and Santos [21]. This implies that the scaling form /(c) in free 
cooling = 0) is identical to the NESS distribution / (c) of the same sys- 
tem, driven by a Gaussian thermostat, provided both forms are rescaled to 
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the same constant width J dcc^f (c) = d/2. This also implies that the scaling 
form for the free case can be measured by performing Monte Carlo simula- 
tions in a steady state [20,21]. The same idea of systematically rescaling the 
velocities has also been used in molecular dynamics simulations of a freely 
cooling system of N inelastic hard spheres [53] . 

Universality: 

The scaling equations (53) for / (c) in systems driven by a Gaussian thermo- 
stat and by white noise, have a universal form, because tJ72 is independent of 
the friction constant 7, the noise strength D and the width vq (00), which may 
depend on the initial distribution. The scaling equations reduce for = 1 to 
those for hard spheres [17], and for v = to those for Maxwell models [35]. 

Perturbative approach: 

The moment of the collision term 1772 = rf7oKiy+2 in the scaling equations 
(53), and the cooling rate C,y = 7oK;y+2^''o in (48) contain for all models the 
quantity Kj^, as given by (38). It depends on the unknown function /(c), 
except for = where K2 = 1 and W2 = c?7o = | (l — a^). For the case of 
inelastic hard spheres, a perturbative method has been developed in [11,17] 
for small inelasticities to solve both integral equations in (53) by expanding 
/ (c) in a series of Sonine polynomials Sp (c^) , i.e. 

00 

f{c) = <P{c){l + Y,apSp{c^)}, (55) 

p=2 

where (c) = tt"''/^ exp (— c^) is the Maxwellian. Here a2 is essentially the 
fourth cumulant which has been calculated explicitly in [17]. For dimensions 
> 1 it turns out to be proportional to the inelasticity (l — o:^). However, 
the one-dimensional case is exceptional because 02 approaches a finite value 
for a ^ 1 [18]. The same method has been successfully applied by Cercignani 
et al. [26] to an inelastic Maxwell model, and by several authors [20-22,50] 
to inelastic hard spheres and related problems. The method can be applied 
to the inelastic IRS-models as well. 

The method focuses on the lower moments of / (c), and the polynomial 
approximation, cut off after 5*2 (c^) , gives a fair representation of / (c) for 
velocities in the thermal range, c < 2. The expansion (55) can then be used 
to calculate Ki, in (52) in the form K,y = k" + (1 — a^) kJ,+- • • . For more details 
we refer to [17]. As an illustration we calculate the lowest approximation 
to in (38) using / (c) ~ (c). The result is, 

K° = y y dcdci Igwl" (j) (c) (t> (ci) 
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The lowest approximation to (54) is then = djoi^^- The damping rate 
in (43) of the r.m.s. velocity becomes then 'yv = (y + 1) 70, in agreement with 
the results for = 1 in [17]. 

4.4 High energy tails 

The polynomial expansion (55) describes / (c) only in the thermal range, but 
contains no meaningful information about velocities c in the asymptotic range 
(c 1). However the high energy tails in the IRS-models with i/ > can be 
determined by a procedure similar to the one used successfully for inelastic 
hard spheres systems [17], as well as in inelastic Maxwell models driven by 
white noise [35] . To do so we make the ansatz of stretched Gaussian behavior 
for the high energy tail, i.e. / (c) ~ B exp [—/3c''] with < 6 < 2 and /? > 0, 
and determine h and (3 by inserting this ansatz in the scaling equation (53) 
and requiring self-consistency. The border line case b ^ 2 corresponds to 
Gaussian tails, and 6 — > suggests power law tails with negative exponents. 

An estimate of the rescaled collision term / (c| /) in (24) is made in [17,35]. 
This suggests that the loss term is asymptotically dominant over the gain 
term as long as the exponent b in exp [— /Sc*"] is restricted to 6 > 0. This 
estimate also applies to the IRS-models as long as the exponent v > and 
the inelasticity 70 is non-vanishing. So the gain term is neglected. Moreover 
the loss term can be simplified in the asymptotic velocity range. Its dominant 
contribution comes from collisions with particles having velocities ci that are 
typically in the thermal range (ci = C (1)). Consequently the collision rate 
|(c — ci) • n\" for asymptotic dynamics may be replaced by c" |c • n\" , and 
the total collision term simplifies to, 

^W/)~AossW/)~-cW(c), (57) 

with given by (40). After these preparations we insert (57) and the stretched 
Gaussian form, /(c) ^ exp [— /3c^] , into the Boltzmann equation (53), and 
follow the procedure sketched above. This gives the following universal results 
(see item 2 below (54)) for the asymptotic high energy tail in d— dimensional 
inelastic models, / (c) exp [—/3c''] with 

6 = inconsistent (GT; u = Q) 

h = v P = dj3^/vw2 (GT; > 0) (53) 

6=1(^ + 2) P=^2\l^ (WN; v>0), 

where the scaling functions and high energy tails of GT-driven and freely 
cooling systems are equivalent. 

We conclude this subsection with some comments. 

driven system,s: 

For < 6 < 2 both GT- and WN-driving lead to consistent asymptotic 
solutions of the scaling equations for 1/ > with overpopulated high energy 
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tails of stretched Gaussian type. We also note, the larger the interactions 

(larger i^-valucs) at large impact energies, the smaller the overpopulation of 
the high energy tails. For v = 1, corresponding to the inelastic version of 
the Very Hard Particle model (see end of Section 3.2) we obtain Gaussian 
behavior (6 = 2) for both types of thermostats, and there arc no longer 
over-populated tails. In the case of stretched Gaussian tails all moments 
J dcc^f (c) < oo. This would not be the case for power law tails. 

In the case of white noise driving, the above results with 6 = 1 + ^v, 
include inelastic hard spheres {v = 1), as well as Maxwell models {v = 0), 
and the results coincide with the detailed predictions for i/ = 1 [17] and 
1^ = [35]. For GT-drivcn or freely cooling IRS-models with v > the tail 
distributions have a stretching exponent b — v. The property j3b ^ I/70 as 
u I Q would suggest power law behavior at the stability threshold, 1/ = 
(Maxwell models), with an exponent a = I/70 (see however next item). The 
result (58) also includes the exponential decay, /(c) ~ exp[— /3c] for inelastic 
hard spheres, found in [17]. 

Stability thresholds: 

In the previous comments we have speculated on power law tails, in agreement 
with the exact solution [30] for the one-dimensional freely cooling Maxwell 
models with / (c) ~ c~^ for c » 1. As will be shown in section 5, the scaling 
form for d-dimcnsional Maxwell models in free cooling does indeed have an 
algebraic tail /(c) ~ l/c''+", but with an exponent a{a) ^ l/7o- However, 
neither the existence of algebraic tails, nor the possible value of the power law 
exponents can be obtained legitimately from the results (58). The reason is 
that the result is not consistent with the a priori assumption /g^jjj ^ -^loss 
at c l^- 1. In fact, gain and loss terms are of the same order of magnitude for 
v = Q. 

Maxwell models as approximations to hard spheres: 

As discussed in item 2 of Section 3.5, inelastic Maxwell models have also been 
introduced as a sensibly looking mathematical simplification of the Boltz- 
mann equation for inelastic hard spheres. The results of the analysis in the 
previous section show that the effect of this simplification on the shape of 
the tail may be very drastic. In the free cooling case the hard sphere tail is 
exponential, / (c) ^ exp [—/3c], whereas in the Maxwell model it is a power 
law tail, /(c) ^ l/c''+". A smaller difference exists in the WN-drivcn case, 
where the hard sphere tail is / (c) ~ exp [— /3c'^/^] , and the Maxwell tail is 
/(c) ~ exp [-/3c]. 

MC- simulations at high and low inelasticity: 

The results (58) predict more than the stretching exponent b. In fact the 
coefhcient /3 for > in / (c) ~ exp [—/3c''] is given in terms of W2 in 
(54). For these z^— values we have described a perturbative calculation for 
tU2, which converges rapidly for a — > 1, but gives poor results for a < 0.6 
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[17,21]. Nevertheless, it is possible to test the predictions of (58) for the GT- 
and the WN-thermostats, inchiding the coefficient /? for all vahics of the 
restitution coefficient a. This can be done by measuring the mean square 
velocity Vq (oo) in MC simulations on driven systems in a NESS for a given 
7 or D. Prom these data 072 can be calculated using (54). 

5 Inelastic Maxwell models 
5.1 Fourier transform method 

As discussed in Section 4.3, the behavior of the high energy tails of the scahng 
form in the IRS-models is only controlled by the loss term in the Boltzmann 
colhsion term, and generates the stretched Gaussian tails. However, in the 
borderline case {1/ = 0) the tail behavior is determined by an interplay be- 
tween gain and loss term, which leads to algebraic tails. This makes the 
analysis more complicated. In this article we only discuss the borderline case 
in free cooling, formed by the Maxwell models, and we refer to [29] for a more 
comprehensive discussion. 

In this section we want to demonstrate that the Boltzmann equation (24) 
for inelastic Maxwell models {p = 0), 

dtF {v, t) = -F{v,t) + ^J J dwF {v**, t) F {w**, t) , (59) 

has a scaling solution with a power law tail. To do so we first consider 
the Fourier transform of the distribution function, <P(fc,t) = (exp [—ik ■ v\), 
which is the characteristic function or generating function of the velocity 
moments. Because F {v, t) is isotropic, ^ (fc, t) is isotropic as well. 

As an auxiliary step we first Fourier transform the gain term in (59), i.e. 

J dv exp [—ik ■ v] /gain {v\F) = 

/ 1 dvdwexp[-ik-v*]F{v,t)F{w,t)= 1 ^{kri+,t)${kr]-,t). (60) 

The transformation needed to obtain the first equality is the same as in (37). 
Then we use (14) to write the exponent as A; • = fe_ • + fe+ • V2, where 

fc+ = pk - nn = kp (jt ■ n^ = kr]^ (n) 

I ^''^ 

fe_ = fe — A;+ [k-[ = k\ 1 — z ik ■ nj = kr]- (n) , 

with p = 1 — q = ^[1 + a) and z = 1 — q^. The Fourier transform of (59) 
then becomes, 

dt$ {k, t) = {k, t)+ f $ {kT,+ (n) , t) $ {kr,- (n) , t) , (62) 

J n 
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where ^ (0, t) = 1 because of (25). In one-dimension r/+ (n) = p and rj- (n) = 
q, and can be replaced by unity. Moreover this equation simplifies to [24] , 

dt^ {k, t)=$ {pk, t) ^ {qk, t) - ^ {k, t) . (63) 

Because F {v, t) is isotropic, only the even moments are non-vanishing, and 
the moment expansion of the characteristic function takes the form, 

^(^'^) = E'^^ ((^ • ^)") = E' (-^^)" ^» w ' (64) 

n ' n 

where the prime indicates that n = even, and the moment m„ {t) is defined 
as, 

TOn {t) = Hn {V^ ln\, (65) 

where /?„ — J^(j^ ■ is given in (40). Moreover, the normalizations (25) 

give Too (t) = 1 and TO2 (i) = ^(32 (v"^) = jv'^- 

Maxwell models have the unusual property that the system of moment 
equations for to„ {t) is closed, and can be solved sequentially. The reason is 
that the rate rh.„ is only a function of lower moments nig (t) with s < n. 

The moment equations are readily obtained by inserting the expansion 
(64) in the Fourier transformed Boltzmann equation (62), and equating the 
coefiicients of equal powers of k. The result is. 



71-2 

mi irin-i (n > 2) , (66) 

1=2 



where all labels {n, I, s} take even values only. Following [33] or Appendix A 
of [36] the functions can be calculated for real positive values of £ and s with 
the result, 

s l+l l+d , 
— ; — ^ z 



h {I, s)=j rj\ (n) Tj'_ (n) = p'(3i 2F1 (^-| 

Xs=l-h {s, 0)-h (0, s)= [ [1 - r?; (n) - r?i (n)] (67) 

=l-p-^P.-.F.{-t,y-\l-q^ 

Here 2-^1(0:, /?; 7 | 2) is a hypergeometric function, and A2 = 2pq/d = 270 is 
easily obtained from the above results. 

Next we consider the Fourier transform of the scaling relation (26), yield- 
ing <P{k,t) = 4'{kv{){t)) with Wo (0 = (0) cxp (— 7ot) according to (41). 
Inserting # in (62) gives the integral equation for the scaling form 0(fc), 
which reads 

-7ofc^<^ {k) + <l>{k) = Jj {krj+) (fcr/_) . (68) 
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For d = 1 it reduces to, 

-pqk^<k{k)+4>{k) = ct>{pk)4>{.qk). (69) 

The scaling form 4'{k) is the generating function for the moments of /(c), 
i.e. 

n (70) 

~1 - + A:4/X4 - fcVe + • • • , 

where n = even, /io = 1 and /i2 = ^/32 (c^) = 1/4 on account of the normal- 
izations (27) and P2 = ^/d (see (40)). By inserting (70) into (68) one obtains 
the recursion relation. 



^ n-2 

l^n = - i-^y^h{l,n-l)niiin-i (n>2) (71) 



with initialization fj,2 = 1/4 and all labels n,l = even. How these moments 
behave as a function of a has been calculated in [36] by numerically solving 
the recursion relation. 



5.2 Small-fe singularity of the characteristic function 

The asymptotic analysis in Section 4.4 of the high energy tail / (c) ~ exp [— /Jc*"] 
for Maxwell models {u = 0) yields 6 = 0, which suggests a power law tail 
/ (c) ~ l/c''+° as the leading large— c behavior in d— dimensional systems. If 
this is indeed the case, then the moments /x„ of the scaling form / (c) are 
convergent if n < o and are divergent if n > a. As we are interested in phys- 
ical solutions which can be normalized, and have a finite energy, a possible 
value of the power law exponent must obey a > 2. 

The characteristic function is in fact a very suitable tool for investigating 
this problem. Suppose the moment /Lt„ with n > a diverges, then the n-th 
derivative of the corresponding generating function also diverges at A; = 0, 
i.e. (f) (k) has a singularity at fc = 0. Then a simple rescaling argument of the 
inverse Fourier transform shows that (j) (k) has a dominant small-fc singularity 
of the form (f> (k) ^ l/cj", where a 7^ even. On the other hand, whenever the 
exponent b in f (c) ^ exp [— /^c**] is positive, - as is the case for the inelastic 
Maxwell model driven by white noise, where 6 = 1 and /3 = ^/2jpq - , then 
all moments are finite, and the characteristic function (j) (fc) is regular at the 
origin, i.e. can be expanded in powers of k^. 

We first illustrate our analysis for the one-dimensional case. As the re- 
quirement of finite total energy imposes the lower bound a > 2 on the ex- 
ponent, we make the ansatz, consistent with (70) that the dominant small-fc 
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singularity has the form, 

= ip+^i^r, (72) 

insert this in (69), and equate the coefficients of equal powers of k. This 
yields, 

iaA2 = Xa or apq=l-p'^- q"-. (73) 

The equation has two roots, a = 2,3, of which a = 3 is the one larger than 
2. Here A is left undetermined. Consequently the one-dimensional scaling 
solution has a power law tail, /(c) ~ 1/c*, in agreement with the exact 
solution in [30]. 

For general dimension we proceed in the same way as in the one- dimensional 
case, insert the ansatz (72) into (68), and equate the coefficients of equal pow- 
ers of k. This yields for the coefficient of fc^ the identity 270 = A2, and for 
the coefficient of k°- the transcendental equation, 

iaA2 =X, = J^[l-nl- r-] ■ (74) 

The equation above obviously has the solution a = 2. We arc however inter- 
ested in the solution with a > 2. In the elastic limit (a ^ 1) the solution is 
simple. There 70 — > and a diverges. The contributions of r]%. on the right 
hand side vanish because r?± < 1, and the result is, 

70 l-a^ 

For general values of a one can conveniently use an integral representation 
of to evaluate Xa and solve the transcendental equation (74) numerically. 
We illustrate the solution method of (74) with the graphical construction in 
Figure 2, where we look for intersections of the line y = \X2 = 70S with the 
curve y = Xg for different values of a. 

The relevant properties of Xg are: (i) limg^o Xg = — 1 ; (ii) Xg is a concave 
function, monotonically increasing with s, and (iii) all eigenvalues for positive 
integers n are positive (see Figure 2). As can be seen from the graphical 
construction, the transcendental equation (74) has two solutions, the trivial 
one (.So = 2) and the solution si = a with a > 2. The numerical solutions for 
d = 2, 3 arc shown in Figure 3 as a function of a, and the Q-dependence of the 
root a (q) can be understood from the graphical construction. In the elastic 
limit an a ] 1 the eigenvalue A2 {a) because of energy conservation. 
In that limit the transcendental equation (74) no longer has a solution with 
a > 2, and a (a) — > 00 according to (75), as it should be. This is consistent 
with a Maxwellian tail distribution in the elastic case. Krapivsky and Ben- 
Naim have in fact solved the transcendental equation asymptotically for large 
d, which gives qualitatively the same results as shown in Figure 3 for two and 
three dimensions. 
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Fig. 2. Graphical solution of Eq.(74) for different values of the parameter a. The 
eigenvalue As is a concave function of s, plotted for different values of the restitution 

coefficient a for the 2-D inelastic Maxwell model. The line y = S70 is plotted for 
a = 0.6,0.8 and a = l(top to bottom). The intersections with As determine the 
points So (filled circles) and si (open circles). Here si = a determines the exponent 
of the power law tail. For the elastic case (a = 1, 70 = 0, energy conservation) there 
is only one intersection point. 




f • ' • ' • ' • ' • ^ 

0.2 0.4 0.6 0.8 1 

a 

Fig. 3. Exponent a (a), which is the non-trivial root of (74), shown as a function 
of the coefficient of restitution a, and denoted in the plot by a. It determines the 
high energy tail 1/c""'"'* of the scaling solution /(c) of the 2-D and 3-D Maxwell 
models. 
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These results establish the existence of scaling solutions /(c) 1/ 0**+" 
with algebraic tails, where the exponent a is the solution of the transcendental 
equation (74) with a > 2. The exponent a (a) behaves as a function of a 
qualitatively the same as in the simple inelastic BGK model of Section 2. 

Using a somewhat different analysis Krapivsky and Ben-Nairn [32] inde- 
pendently obtained the same results for the algebraic tails in freely cooling 
Maxwell models. 



5.3 Beyond asymptotic analysis 

Free cooling: 

As mentioned already in Section 3.3, Baldassarri et al. [30] have obtained an 
exact solution / (c) of the scaling equation for the one-dimensional Maxwell 
model in free cooling, and they demonstrated the importance of this scaling 
solution by means of Monte Carlo simulations. In doing so they found that 
F {v,t) for different classes of initial distributions can be collapsed for long 
times t on this exact scaling solution / (c), when vo (t) F {v,t) is plotted versus 
c = v/vo {t). 

We briefly illustrate here how this solution is obtained from the one- 
dimensional scaling equation (69). One first verifies that the following func- 
tion, 

(/)(fc) = (l + i?|fc|)exp(-i?|fc|), (76) 

with arbitrary positive i? is a solution, that can be Fourier inverted. It gives 
the scaling form, 

1 

TTI? (1 + ^/^^f 

To determine the scaling form that satisfies the normalizations (27), we ex- 
pand (76) in powers of to obtain. 



(A(fc) = l-il?2fc2^1^3|fc|^ + ... , (78) 

Comparison of this result with (70) shows that -d = 1/^/2. Moreover, it 
confirms the ansatz (72) used to find solutions with small-fc singularities. 
Comparison also shows the value of the coefficient A= 1/ [6\/2] in (72), and 
the high energy tail is /(c) ~ A/c with A = 1/ \tt^/2\. The coefficients A 
and A can not be determined within the asymptotic method. 

White noise driving: 

For the one-dimensional Maxwell model driven by white noise the steady state 
solution has also been found exactly with the help of the Fourier transform 
method [24]. The characteristic function satisfies in that case a nonlinear 
finite difference equation, which can be solved by iteration. However, the 
analytic structure is rather complex, and makes it difficult to extract analytic 
information from that solution. 
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The observation that / (c) exp [— /3|c|] has an exponential high energy 
tail, was probably first made in numerical work of van dcr Hart and Nicnhuis 
[38] . A more detailed analytic prediction about the asymptotic tail was given 
in [35], where it was shown that /3 = 1/8/ (1 — a^) for all Maxwell models, 
independent of the dimensionality (d = 1, 2, • • • ), at least with the normaliza- 
tion, (c^) = ^d, used in this article. This result follows directly from (58) and 
the value = djo = 3 (l ~ o:^), given in item 3 below (54). Furthermore, 
additional numerical and analytical work was also published by Marconi and 
Puglisi [54], and by Antal et al. [39]. Only recently more detailed analytic 
results have been extracted from the rather complex structure of the exact 
solution [48,40]. 

6 Conclusions and perspectives 

We have studied asymptotic properties of scaling or similarity solutions, 
F {v, t) = {vo (i))"'' / {v/vo (t)), of the nonlinear Boltzmann equation in spa- 
tially homogeneous systems composed of particles with inelastic interactions 
for large times and large velocities. The large t— and w— scales arc relevant 
because on such scales the universal features of the solutions survive, while 
details of the initial distributions, of interaction strength, and degree of in- 
elasticity, arc; mostly lost. The behavior of these scaling states, which describe 
nonequilibrium steady states (NESS), is less universal than the state of ther- 
mal equilibrium, because the form of the NESS distribution /(c) depends on 
the way of driving the dissipative systems. Scaling solutions are very well 
suited to expose the universal features of the velocity distribution functions, 
because the velocities, c = v/vo{t), are measured in units of the r.m.s. veloc- 
ity or instantaneous width vq (t) of the distribution. 

The real importance of the scaling solutions is that the actual solutions 
F{v, t) for large classes of initial distributions F{v, 0) (essentially initial data 
without over-populated tails -sec [36]) rapidly approach these scaling solu- 
tions in the sense that after a short transient time the data Vq {t) F {v,t) 
can be collapsed on a single scaling form / (c), as first observed by Baldas- 
sarri et al. in MC simulations of the nonlinear Boltzmann equation for a 
one-dimensional inelastic Maxwell model. This rapid approach to a universal 
scaling function applies both to systems driven by a Gaussian or by a white 
noise thermostat, as well as for freely cooling systems, which show scaling 
behavior, identical to systems driven by a Gaussian thermostat. 

Originally this scenario had the status of a conjecture for systems of inelas- 
tic hard spheres (i> = 1) and inelastic Maxwell models(i^ = 0), as formulated 
in [36] where also some analytical evidence for the approach to a scaling form 
has been presented. For Maxwell models the conjecture has been rigorously 
proven in the mean time by Bobylev et al. [37]. The analysis in Section 4.1 
also suggests what the basic criterion us for the approach of IRS-models to 
the scaling form. If the energy balance equation (48) has a stable/attractive 
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fixed point solution, t;o(oo), then the distribution function F{v, t) approaches 

- in the sense detailed in the paragraphs above - a stable NESS described 
by the scaling solution. 

The analysis in Section 4.1 of the energy balance equation also suggests 

that cases of marginal stability are candidates for power law tails, where the 
freely evolving inelastic Maxwell model is a well known example. A further 
candidate system for power law tails is an IRS-model with u = —2, driven 
by white noise. These properties, and the relation between marginal stability 
and power law tails is further explored in [29] , where in addition the analytic 
predictions of the present paper are confirmed by MC simulations. 

In this article we have focussed on the properties of the scaling solu- 
tion, and in particular on its high energy tail. Here we have introduced the 
Boltzmann equation for new classes of inelastic interactions, named Inelastic 
Repulsive Scatterers or IRS-models, corresponding to pseudo-repulsive power 
law potentials, with collision probabilities proportional to A ~ \g ■ n\'^ , cover- 
ing hard scatterers like inelastic hard spheres (v = 1) and soft scatterers like 
pseudo-Maxwell molecules = 0). The energy loss in an inelastic interaction 
is proportional to the inelasticity, 70 ~ (l — a^), where a is the coefficient 
of restitution. We have studied two typical cases: freely evolving systems in 
homogeneous cooling states, without energy supply, and systems with energy 
supply, driven by Gaussian thermostats with negative friction or driven by 
Gaussian white noise. 

The homogeneous cooling laws in these systems are described by the gran- 
ular temperature T = Vq, with a long time decay as i"^/*^ or exp[— 27yT], 
where t is the external (laboratory) time, and r is the internal time or collision 
counter, and 7,^ is proportional to the inelasticity 70. 

An interesting feature of all IRS-models is that their scaling forms gener- 
ically have a stretched Gaussian high energy tail, /(c) ^ exp[— /3c*]. In 
freely cooling IRS-models with v > the stretching exponent h = v satis- 
fies Q <b = V <2. Among the IRS-models at the stability threshold only the 
freely cooling Maxwell model has been analyzed in the present paper, and it 
does yield power law tails, / (c) ~ l/c''+", and the exponent 0, has been cal- 
culated from a transcendental equation. In the IRS-models with v > Q, driven 
by white noise, there exist again stretched Gaussian tails with 6 = ^{v + 2). 
All known results for higher dimensional inelastic hard spheres {v — 0) and 
inelastic Maxwell models {v = 1) are recovered in the present analysis. 

In IRS-models with positive v the tails are always stretched Gaussians, 
and they are determined only by the loss term in the Boltzmann equation. In 
the freely cooling Maxwell model, and more generally at a stability threshold 

- which are the cases leading to power law tails [29] - the loss and gain term 
in the nonlinear Boltzmann equation are of comparable size, and partially 
balancing each other. 

We also have analyzed in Section 2 an inelastic BGK model, introduced by 
Brey et al [44] , and generalized it to an energy dependent collision frequency. 
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and we added energy source terms to this kinetic equation as well. In all cases 
the scaling solution can be calculated exactly. The freely cooling model shows 
an algebraic tail, / (c) ~ l/c''+" with an exponent a ~ I/70 ~ 2/ (l — a^). 
For white noise driving one finds asymptotically / (c) ~ exp [—(3c] with /3 

It is interesting to observe that the overpopulated high energy tails in 
the inelastic BGK-models, both for free cooling as well as for white noise 
driving, arc essentially the same as for the more complicated Maxwell models. 
However, the present extension to an energy dependent collision frequency 
does not capture the generic features of the less-singular stretched Gaussian 
tails for the IRS-models with v > , as predicted by the nonlinear Boltzmann 
equation for these models. 
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